home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / symb_lib / prisa.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-12-11  |  19.3 KB  |  527 lines
  1. /******************************************************************************
  2. * Prisa.c - piecewise ruled srf approximation and layout (prisa).          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Apr. 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include "symb_loc.h"
  8.  
  9. static CagdSrfStruct *SymbPiecewiseRuledSrfAux(CagdSrfStruct *Srf,
  10.                            CagdBType ConsistentDir,
  11.                            CagdRType Epsilon,
  12.                            CagdSrfDirType Dir);
  13. static void SymbInterCircCirc(CagdRType Center1[3],
  14.                   CagdRType Radius1,
  15.                   CagdRType Center2[3],
  16.                   CagdRType Radius2,
  17.                   CagdRType Inter1[3],
  18.                   CagdRType Inter2[3]);
  19.  
  20. /*****************************************************************************
  21. * DESCRIPTION:                                                               M
  22. * Computes a piecewise ruled surface approximation to a given set of         M
  23. * surfaces with given Epsilon, and lay them out "nicely" onto the XY plane,  M
  24. * by approximating each ruled surface as a developable surface with          M
  25. * SamplesPerCurve samples.                             M
  26. *   Dir controls the direction of ruled approximation, SpaceScale and        M
  27. * Offset controls the placement of the different planar pieces.             M
  28. *   Prisa is the hebrew word for the process of flattening out a three       M
  29. * dimensional surface. I have still to find an english word for it.         M
  30. *                                                                            *
  31. * PARAMETERS:                                                                M
  32. *   Srfs:            To approximate and faltten out.                         M
  33. *   SamplesPerCurve: During the approximation of a ruled surface as a        M
  34. *             developable surface.                     M
  35. *   Epsilon:         Accuracy control for the piecewise ruled surface        M
  36. *             approximation.                         M
  37. *   Dir:             Direction of ruled/developable surface approximation.   M
  38. *             Either U or V.                         M
  39. *   Space:           A vector in the XY plane to denote the amount of        M
  40. *             translation from one flattened out surface to the next. M
  41. * RETURN VALUE:                                                              M
  42. *   CagdSrfStruct *:   A list of planar surfaces denoting the layout (prisa) M
  43. *               of the given Srfs to the accuracy requested.         M
  44. *                                                                            *
  45. * KEYWORDS:                                                                  M
  46. *   SymbAllPrisaSrfs, layout, prisa                                          M
  47. *****************************************************************************/
  48. CagdSrfStruct *SymbAllPrisaSrfs(CagdSrfStruct *Srfs,
  49.                 int SamplesPerCurve,
  50.                     CagdRType Epsilon,
  51.                 CagdSrfDirType Dir,
  52.                     CagdVType Space)
  53. {
  54.     int SrfIndex = 1;
  55.     CagdSrfStruct *Srf,
  56.     *PrisaSrfsAll = NULL;
  57.     CagdVType Offset;
  58.  
  59.     for (Srf = Srfs; Srf != NULL; Srf = Srf -> Pnext, SrfIndex++) {
  60.     CagdSrfStruct *TSrf, *RSrf, *RuledSrfs;
  61.  
  62.     if (Epsilon > 0) {
  63.         RuledSrfs = SymbPiecewiseRuledSrfApprox(Srf, FALSE, Epsilon, Dir);
  64.  
  65.         Offset[0] = SrfIndex * Space[0];
  66.         Offset[1] = 0.0;
  67.         Offset[2] = 0.0;
  68.  
  69.         for (RSrf = RuledSrfs; RSrf != NULL; RSrf = RSrf -> Pnext) {
  70.         TSrf = SymbPrisaRuledSrf(RSrf, SamplesPerCurve,
  71.                      Space[1], Offset);
  72.         TSrf -> Pnext = PrisaSrfsAll;
  73.         PrisaSrfsAll = TSrf;
  74.         }
  75.  
  76.         CagdSrfFreeList(RuledSrfs);
  77.     }
  78.     else {
  79.         /* Return the 3D ruled surface approximation. */
  80.         RuledSrfs = SymbPiecewiseRuledSrfApprox(Srf, FALSE, -Epsilon, Dir);
  81.  
  82.         for (TSrf = RuledSrfs;
  83.          TSrf -> Pnext != NULL;
  84.          TSrf = TSrf -> Pnext);
  85.         TSrf -> Pnext = PrisaSrfsAll;
  86.         PrisaSrfsAll = RuledSrfs;
  87.     }
  88.     }
  89.  
  90.     return PrisaSrfsAll;
  91. }
  92.  
  93. /*****************************************************************************
  94. * DESCRIPTION:                                                               M
  95. * Constructs a piecewise ruled surface approximation to the given surface,   M
  96. * Srf, in the given direction, Dir, that is close to the surface to within   M
  97. * Epsilon.                                     M
  98. *   If ConsitentDir then ruled surface parametrization is set to be the      M
  99. * same as original surface Srf. Otherwise, ruling dir is always             M
  100. * CAGD_CONST_V_DIR.                                 M
  101. *   Surface is assumed to have point types E3 or P3 only.             M
  102. *                                                                            *
  103. * PARAMETERS:                                                                M
  104. *   Srf:           To approximate using piecewise ruled surfaces.            M
  105. *   ConsistentDir: Do we want parametrization to be the same as Srf?         M
  106. *   Epsilon:       Accuracy of piecewise ruled surface approximation.        M
  107. *   Dir:           Direction of piecewise ruled surface approximation.       M
  108. *           Either U or V.                         M
  109. *                                                                            *
  110. * RETURN VALUE:                                                              M
  111. *   CagdSrfStruct *:   A list of ruled surfaces approximating Srf to within  M
  112. *                      Epsilon in direction Dir.                             M
  113. *                                                                            *
  114. * KEYWORDS:                                                                  M
  115. *   SymbPiecewiseRuledSrfApprox, layout, prisa, ruled surface approximation  M
  116. *****************************************************************************/
  117. CagdSrfStruct *SymbPiecewiseRuledSrfApprox(CagdSrfStruct *Srf,
  118.                        CagdBType ConsistentDir,
  119.                        CagdRType Epsilon,
  120.                        CagdSrfDirType Dir)
  121. {
  122.     CagdSrfStruct *RuledSrfs;
  123.  
  124.     switch (Srf -> PType) {
  125.     case CAGD_PT_E3_TYPE:
  126.     case CAGD_PT_P3_TYPE:
  127.         Srf = CagdSrfCopy(Srf);
  128.         break;
  129.     case CAGD_PT_P1_TYPE:
  130.     case CAGD_PT_P2_TYPE:
  131.     case CAGD_PT_P4_TYPE:
  132.     case CAGD_PT_P5_TYPE:
  133.         Srf = CagdCoerceSrfTo(Srf, CAGD_PT_P3_TYPE);
  134.         break;
  135.     case CAGD_PT_E1_TYPE:
  136.     case CAGD_PT_E2_TYPE:
  137.     case CAGD_PT_E4_TYPE:
  138.     case CAGD_PT_E5_TYPE:
  139.         Srf = CagdCoerceSrfTo(Srf, CAGD_PT_E3_TYPE);
  140.         break;
  141.     default:
  142.         SYMB_FATAL_ERROR(SYMB_ERR_UNSUPPORT_PT);
  143.         Srf = NULL;
  144.         break;
  145.     }
  146.  
  147.     RuledSrfs = SymbPiecewiseRuledSrfAux(Srf, ConsistentDir, Epsilon, Dir);
  148.     CagdSrfFree(Srf);
  149.     return RuledSrfs;
  150. }
  151.  
  152. /*****************************************************************************
  153. * DESCRIPTION:                                                               *
  154. * Auxiliary function to function SymbPiecewiseRuledSrfApprox             *
  155. *****************************************************************************/
  156. static CagdSrfStruct *SymbPiecewiseRuledSrfAux(CagdSrfStruct *Srf,
  157.                            CagdBType ConsistentDir,
  158.                            CagdRType Epsilon,
  159.                            CagdSrfDirType Dir)
  160. {
  161.     CagdSrfStruct *DiffSrf, *DistSqrSrf,
  162.     *RuledSrf = CagdSrfCopy(Srf);
  163.     CagdRType *XPts, *WPts, UMin, UMax, VMin, VMax,
  164.     **Points = RuledSrf -> Points,
  165.     MaxError = 0.0,
  166.     t = 0.0;
  167.     int i, j, k, Index,
  168.     ULength = RuledSrf -> ULength,
  169.     VLength = RuledSrf -> VLength;
  170.     PointType E3PtStart, E3PtEnd, E3Pt;
  171.  
  172.     switch (Dir) {
  173.     case CAGD_CONST_U_DIR:
  174.         for (j = 0; j < VLength; j++) {
  175.         /* First order approximation to the ratios of the   */
  176.         /* distance of interior point to the end points.    */
  177.         CagdCoerceToE3(E3PtStart, Points,
  178.                    CAGD_MESH_UV(RuledSrf, ULength / 2, 0),
  179.                    Srf -> PType);
  180.         CagdCoerceToE3(E3PtEnd, Points,
  181.                    CAGD_MESH_UV(RuledSrf, ULength / 2,
  182.                                   VLength - 1),
  183.                    Srf -> PType);
  184.         CagdCoerceToE3(E3Pt, Points,
  185.                    CAGD_MESH_UV(RuledSrf, ULength / 2, j),
  186.                    Srf -> PType);
  187.         PT_SUB(E3PtStart, E3PtStart, E3PtEnd);
  188.         PT_SUB(E3Pt, E3Pt, E3PtEnd);
  189.         t = PT_LENGTH(E3PtStart);
  190.         t = APX_EQ(t, 0.0) ? 0.5 : PT_LENGTH(E3Pt) / t;
  191.  
  192.         for (i = 0; i < ULength; i++) {
  193.             CagdCoerceToE3(E3PtStart, Points,
  194.                    CAGD_MESH_UV(RuledSrf, i, 0),
  195.                    Srf -> PType);
  196.             CagdCoerceToE3(E3PtEnd, Points,
  197.                    CAGD_MESH_UV(RuledSrf, i, VLength - 1),
  198.                    Srf -> PType);
  199.             PT_BLEND(E3Pt, E3PtStart, E3PtEnd, t);
  200.  
  201.             Index = CAGD_MESH_UV(RuledSrf, i, j);
  202.             if (CAGD_IS_RATIONAL_PT(RuledSrf -> PType)) {
  203.             for (k = 0; k < 3; k++)
  204.                 Points[k + 1][Index] = 
  205.                 E3Pt[k] * Points[0][Index];
  206.             }
  207.             else {
  208.             for (k = 0; k < 3; k++)
  209.                 Points[k + 1][Index] = E3Pt[k];
  210.             }
  211.         }
  212.         }
  213.         break;
  214.     case CAGD_CONST_V_DIR:
  215.         for (i = 0; i < ULength; i++) {
  216.         /* First order approximation to the ratios of the   */
  217.         /* distance of interior point to the end points.    */
  218.         CagdCoerceToE3(E3PtStart, Points,
  219.                    CAGD_MESH_UV(RuledSrf, 0, VLength / 2),
  220.                    Srf -> PType);
  221.         CagdCoerceToE3(E3PtEnd, Points,
  222.                    CAGD_MESH_UV(RuledSrf, ULength - 1,
  223.                               VLength / 2),
  224.                    Srf -> PType);
  225.         CagdCoerceToE3(E3Pt, Points,
  226.                    CAGD_MESH_UV(RuledSrf, i, VLength / 2),
  227.                    Srf -> PType);
  228.         PT_SUB(E3PtStart, E3PtStart, E3PtEnd);
  229.         PT_SUB(E3Pt, E3Pt, E3PtEnd);
  230.         t = PT_LENGTH(E3PtStart);
  231.         t = APX_EQ(t, 0.0) ? 0.5 : PT_LENGTH(E3Pt) / t;
  232.  
  233.         for (j = 0; j < VLength; j++) {
  234.             CagdCoerceToE3(E3PtStart, Points,
  235.                    CAGD_MESH_UV(RuledSrf, 0, j),
  236.                    Srf -> PType);
  237.             CagdCoerceToE3(E3PtEnd, Points,
  238.                    CAGD_MESH_UV(RuledSrf, ULength - 1, j),
  239.                    Srf -> PType);
  240.             PT_BLEND(E3Pt, E3PtStart, E3PtEnd, t);
  241.  
  242.             Index = CAGD_MESH_UV(RuledSrf, i, j);
  243.             if (CAGD_IS_RATIONAL_PT(RuledSrf -> PType)) {
  244.             for (k = 0; k < 3; k++)
  245.                 Points[k + 1][Index] = 
  246.                 E3Pt[k] * Points[0][Index];
  247.             }
  248.             else {
  249.             for (k = 0; k < 3; k++)
  250.                 Points[k + 1][Index] = E3Pt[k];
  251.             }
  252.         }
  253.         }
  254.         break;
  255.     default:
  256.         SYMB_FATAL_ERROR(SYMB_ERR_DIR_NOT_CONST_UV);
  257.         break;
  258.     }
  259.  
  260.     DiffSrf = SymbSrfSub(Srf, RuledSrf);
  261.     CagdSrfFree(RuledSrf);
  262.     DistSqrSrf = SymbSrfDotProd(DiffSrf, DiffSrf);
  263.     CagdSrfFree(DiffSrf);
  264.     XPts = DistSqrSrf -> Points[1];
  265.     WPts = CAGD_IS_RATIONAL_PT(DistSqrSrf -> PType) ?
  266.                     DistSqrSrf -> Points[0] : NULL;
  267.  
  268.     for (i = DistSqrSrf -> ULength * DistSqrSrf -> VLength; i > 0; i--) {
  269.     CagdRType
  270.         V = WPts != NULL ? *XPts++ / *WPts++ : *XPts++;
  271.  
  272.     if (MaxError < V)
  273.         MaxError = V;
  274.     }
  275.     CagdSrfFree(DistSqrSrf);
  276.  
  277.     if (MaxError > SQR(Epsilon)) {
  278.     /* Subdivide and try again. */
  279.     CagdSrfStruct *Srf1, *Srf2, *RuledSrf1, *RuledSrf2;
  280.  
  281.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  282.     t = Dir == CAGD_CONST_V_DIR ? (UMax + UMin) / 2 : (VMax + VMin) / 2;
  283.     Srf1 = CagdSrfSubdivAtParam(Srf, t, CAGD_OTHER_DIR(Dir));
  284.     Srf2 = Srf1 -> Pnext;
  285.     Srf1 -> Pnext = NULL;
  286.  
  287.     RuledSrf1 = SymbPiecewiseRuledSrfAux(Srf1, ConsistentDir, Epsilon, Dir);
  288.     RuledSrf2 = SymbPiecewiseRuledSrfAux(Srf2, ConsistentDir, Epsilon, Dir);
  289.     CagdSrfFree(Srf1);
  290.     CagdSrfFree(Srf2);
  291.  
  292.     for (Srf1 = RuledSrf1; Srf1 -> Pnext != NULL; Srf1 = Srf1 -> Pnext);
  293.     Srf1 -> Pnext = RuledSrf2;
  294.     return RuledSrf1;
  295.     }
  296.     else {
  297.     /* Returns the ruled surface as a linear in the ruled direction. */
  298.     CagdCrvStruct
  299.         *Crv1 = CagdCrvFromMesh(Srf, 0, CAGD_OTHER_DIR(Dir)),
  300.         *Crv2 = CagdCrvFromMesh(Srf,
  301.                     Dir == CAGD_CONST_V_DIR ? ULength - 1
  302.                                 : VLength - 1,
  303.                     CAGD_OTHER_DIR(Dir));
  304.  
  305.     RuledSrf = CagdRuledSrf(Crv1, Crv2, 2, 2);
  306.     if (ConsistentDir && Dir == CAGD_CONST_V_DIR) {
  307.         /* Needs to reverse the ruled surface so it matches Srf. */
  308.         CagdSrfStruct
  309.         *TSrf = CagdSrfReverse2(RuledSrf);
  310.  
  311.         CagdSrfFree(RuledSrf);
  312.         RuledSrf = TSrf;
  313.     }
  314.  
  315.     if (CAGD_IS_BSPLINE_SRF(Srf)) {
  316.         /* Updates the knot vector domain. */
  317.         CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  318.         if (Dir == CAGD_CONST_V_DIR)
  319.             BspKnotAffineTrans(RuledSrf -> UKnotVector,
  320.                    RuledSrf -> ULength + RuledSrf -> UOrder,
  321.                    UMin, UMax - UMin);
  322.         else
  323.             BspKnotAffineTrans(RuledSrf -> VKnotVector,
  324.                    RuledSrf -> VLength + RuledSrf -> VOrder,
  325.                    VMin, VMax - VMin);
  326.     }
  327.     CagdCrvFree(Crv1);
  328.     CagdCrvFree(Crv2);
  329.     return RuledSrf;
  330.     }
  331. }
  332.  
  333. /*****************************************************************************
  334. * DESCRIPTION:                                                               M
  335. * Layout a single ruled surface, by approximating it as a set of polygons.   M
  336. *   The given ruled surface might be non-developable, in which case         M
  337. * approximation will be of a surface with no twist.                 M
  338. *   The ruled surface is assumed to be constructed using CagdRuledSrf and    M
  339. * that the ruled direction is consisnt and is always CAGD_CONST_V_DIR.         M
  340. *                                                                            *
  341. * PARAMETERS:                                                                M
  342. *   Srf:             A ruled surface to layout flat on the XY plane.         M
  343. *   SamplesPerCurve: During the approximation of a ruled surface as a        M
  344. *             developable surface.                     M
  345. *   Space:         Increment on Y on the offset vector, after this         M
  346. *             surface was placed in the XY plane.             M
  347. *   Offset:         A vector in the XY plane to denote the amount of        M
  348. *             translation for the flatten surface in the XY plane.    M
  349. *                                                                            *
  350. * RETURN VALUE:                                                              M
  351. *   CagdSrfStruct *:  A planar surface in the XY plane approximating the     M
  352. *                     falttening process of Srf.                 M
  353. *                                                                            *
  354. * KEYWORDS:                                                                  M
  355. *   SymbPrisaRuledSrf, layout, prisa                                         M
  356. *****************************************************************************/
  357. CagdSrfStruct *SymbPrisaRuledSrf(CagdSrfStruct *Srf,
  358.                  int SamplesPerCurve,
  359.                  CagdRType Space,
  360.                  CagdVType Offset)
  361. {
  362.     int i, j,
  363.     VLength = Srf -> VLength;
  364.     CagdRType Angle, PtTmp1[3], PtTmp2[3], PtTmp3[3];
  365.     CagdCrvStruct
  366.     *Crv1 = CagdCrvFromMesh(Srf, 0, CAGD_CONST_V_DIR),
  367.     *Crv2 = CagdCrvFromMesh(Srf, VLength - 1, CAGD_CONST_V_DIR);
  368.     CagdPolylineStruct
  369.     *Poly1 = SymbCrv2Polyline(Crv1, SamplesPerCurve, FALSE, TRUE),
  370.     *Poly2 = SymbCrv2Polyline(Crv2, SamplesPerCurve, FALSE, TRUE),
  371.     *Poly1Prisa = CagdPolylineNew(Poly1 -> Length),
  372.     *Poly2Prisa = CagdPolylineNew(Poly2 -> Length);
  373.     CagdPolylnStruct
  374.         *Pt1 = Poly1 -> Polyline,
  375.         *Pt2 = Poly2 -> Polyline,
  376.         *Pt1Prisa = Poly1Prisa -> Polyline,
  377.         *Pt2Prisa = Poly2Prisa -> Polyline;
  378.     CagdMType Mat1, Mat;
  379.     CagdBBoxStruct BBox;
  380.  
  381.     CagdCrvFree(Crv1);
  382.     CagdCrvFree(Crv2);
  383.  
  384.     /* Anchor the location of the first line. */
  385.     for (j = 0; j < 3; j++) {
  386.     Pt1Prisa -> Pt[j] = 0.0;
  387.     Pt2Prisa -> Pt[j] = 0.0;
  388.     }
  389.     PT_SUB(PtTmp1, Pt1 -> Pt, Pt2 -> Pt);
  390.     Pt2Prisa -> Pt[1] = PT_LENGTH(PtTmp1);
  391.  
  392.     /* Move alternately along Poly1 and Poly2 and anchor the next point of */
  393.     /* the next triangle using the distances to current Pt1 and Pt2.       */
  394.     for (i = 2; i < Poly1 -> Length + Poly2 -> Length; i++) {
  395.     CagdRType Dist1, Dist2, Inter1[3], Inter2[3],
  396.         *NextPt = i & 0x01 ? Pt1[1].Pt : Pt2[1].Pt;
  397.  
  398.     /* Compute distance from previous two locations to new location. */ 
  399.     PT_SUB(PtTmp1, Pt1 -> Pt, NextPt);
  400.     Dist1 = PT_LENGTH(PtTmp1);
  401.     PT_SUB(PtTmp1, Pt2 -> Pt, NextPt);
  402.     Dist2 = PT_LENGTH(PtTmp1);
  403.  
  404.     /* Find the (two) intersection points of circles with radii Dist?. */
  405.     SymbInterCircCirc(Pt1Prisa -> Pt, Dist1, Pt2Prisa -> Pt, Dist2,
  406.               Inter1, Inter2);
  407.  
  408.     /* Find which of the two intersection points is the "good" one. */
  409.     PT_SUB(PtTmp1, Inter1, Pt1Prisa -> Pt);
  410.     PT_SUB(PtTmp2, Inter1, Pt2Prisa -> Pt);
  411.     CROSS_PROD(PtTmp3, PtTmp1, PtTmp2);
  412.     if (i & 0x01) {
  413.         Pt1Prisa++;
  414.         for (j = 0; j < 3; j++)
  415.         Pt1Prisa -> Pt[j] = (PtTmp3[2] > 0 ? Inter2[j] : Inter1[j]);
  416.         Pt1++;
  417.     }
  418.     else {
  419.         Pt2Prisa++;
  420.         for (j = 0; j < 3; j++)
  421.         Pt2Prisa -> Pt[j] = (PtTmp3[2] > 0 ? Inter2[j] : Inter1[j]);
  422.         Pt2++;
  423.     }
  424.     }
  425.  
  426.     /* Save centering location so we can orient the resulting data nicely. */
  427.     PT_COPY(PtTmp1, Poly1Prisa -> Polyline[Poly1Prisa -> Length / 2].Pt);
  428.     PT_COPY(PtTmp2, Poly2Prisa -> Polyline[Poly2Prisa -> Length / 2].Pt);
  429.     PT_SUB(PtTmp3, PtTmp2, PtTmp1);
  430.  
  431.     Crv1 = CnvrtPolyline2LinBsplineCrv(Poly1Prisa);
  432.     Crv2 = CnvrtPolyline2LinBsplineCrv(Poly2Prisa);
  433.     CagdPolylineFree(Poly1);
  434.     CagdPolylineFree(Poly2);
  435.     CagdPolylineFree(Poly1Prisa);
  436.     CagdPolylineFree(Poly2Prisa);
  437.  
  438.     Srf = CagdRuledSrf(Crv1, Crv2, 2, 2);
  439.     CagdCrvFree(Crv1);
  440.     CagdCrvFree(Crv2);
  441.  
  442.     /* Translate PtTmp1 to the origin. */
  443.     MatGenMatTrans(-PtTmp1[0], -PtTmp1[1], -PtTmp1[2], Mat);
  444.  
  445.     /* Rotate PtTmp3 direction to the +Y direction. */
  446.     Angle = atan2(PtTmp3[1], PtTmp3[0]);
  447.     MatGenMatRotZ1(M_PI / 2 - Angle, Mat1);
  448.  
  449.     MatMultTwo4by4(Mat, Mat, Mat1);
  450.  
  451.     CagdSrfMatTransform(Srf, Mat);
  452.  
  453.     /* Translate by the Offset. */
  454.     CagdSrfBBox(Srf, &BBox);
  455.     MatGenMatTrans(Offset[0], Offset[1] - BBox.Min[1], Offset[2], Mat);
  456.     Offset[1] += (BBox.Max[1] - BBox.Min[1]) + Space;
  457.  
  458.     CagdSrfMatTransform(Srf, Mat);
  459.  
  460.     return Srf;
  461. }
  462.  
  463. /*****************************************************************************
  464. * DESCRIPTION:                                                               *
  465. * Finds the two intersection points of the given two planar circles. It is   *
  466. * assumed intersection exists.                             *
  467. *                                                                            *
  468. * PARAMETERS:                                                                *
  469. *   Center1, Radius1:  Geometry of first circle.                             *
  470. *   Center2, Radius2:  Geometry of second circle.                            *
  471. *   Inter1, Inter2:    Where the two intersection locations will be placed.  *
  472. *                                                                            *
  473. * RETURN VALUE:                                                              *
  474. *   void                                                                     *
  475. *****************************************************************************/
  476. static void SymbInterCircCirc(CagdRType Center1[3],
  477.                   CagdRType Radius1,
  478.                   CagdRType Center2[3],
  479.                   CagdRType Radius2,
  480.                   CagdRType Inter1[3],
  481.                   CagdRType Inter2[3])
  482. {
  483.     CagdRType A, B, C, Delta,
  484.     a = Center2[0] - Center1[0],
  485.     b = Center2[1] - Center1[1],
  486.     c = (SQR(Radius1) - SQR(Radius2) +
  487.          SQR(Center2[0]) - SQR(Center1[0]) +
  488.          SQR(Center2[1]) - SQR(Center1[1])) / 2;
  489.  
  490.     if (PT_APX_EQ(Center1, Center2)) {
  491.     Inter1[0] = Inter2[0] = Radius1;
  492.     Inter1[1] = Inter2[1] = 0.0;
  493.     }
  494.     else if (ABS(a) > ABS(b)) {
  495.     /* Solve for Y first. */
  496.     A = SQR(b) / SQR(a) + 1;
  497.     B = 2 * (b * Center1[0] / a - b * c / SQR(a) - Center1[1]);
  498.     C = SQR(c / a) + SQR(Center1[0]) + SQR(Center1[1])
  499.                 -2 * c * Center1[0] / a - SQR(Radius1);
  500.     Delta = SQR(B) - 4 * A * C;
  501.     if (Delta < 0)  /* If no solution, do something almost reasonable. */
  502.         Delta = 0;
  503.     Inter1[1] = (-B + sqrt(Delta)) / (2 * A);
  504.     Inter2[1] = (-B - sqrt(Delta)) / (2 * A);
  505.  
  506.     Inter1[0] = (c - b * Inter1[1]) / a;
  507.     Inter2[0] = (c - b * Inter2[1]) / a;
  508.     }
  509.     else {
  510.     /* Solve for X first. */
  511.     A = SQR(a) / SQR(b) + 1;
  512.     B = 2 * (a * Center1[1] / b - a * c / SQR(b) - Center1[0]);
  513.     C = SQR(c / b) + SQR(Center1[0]) + SQR(Center1[1])
  514.                 -2 * c * Center1[1] / b - SQR(Radius1);
  515.     Delta = SQR(B) - 4 * A * C;
  516.     if (Delta < 0)  /* If no solution, do something almost reasonable. */
  517.         Delta = 0;
  518.     Inter1[0] = (-B + sqrt(Delta)) / (2 * A);
  519.     Inter2[0] = (-B - sqrt(Delta)) / (2 * A);
  520.  
  521.     Inter1[1] = (c - a * Inter1[0]) / b;
  522.     Inter2[1] = (c - a * Inter2[0]) / b;
  523.     }
  524.  
  525.     Inter1[2] = Inter2[2] = 0.0;
  526. }
  527.